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We analyse the efficiency of several simulation methods which we have recently proposed for 
calculating rate constants for rare events in stochastic dynamical systems, in or out of equilibrium. 
We derive analytical expressions for the computational cost of using these methods, and for the 
statistical error in the final estimate of the rate constant, for a given computational cost. These 
expressions can be used to determine which method to use for a given problem, to optimize the choice 
of parameters, and to evaluate the significance of the results obtained. We apply the expressions 
to the two-dimensional non-equilibrium rare event problem proposed by Maier and Stein. For this 
problem, our analysis gives accurate quantitative predictions for the computational efficiency of the 
three methods. 



I. INTRODUCTION 

Rare events are processes which happen rapidly, yet in- 
frequently. Specialized techniques are required in order 
to study these events using computer simulation. This 
is because, in "brute force" simulations, the vast ma- 
jority of the computational effort is used in simulating 
the uninteresting waiting periods between events, so that 
observing enough events for reliable statistical analysis 
is generally impossible. The quantities of interest from 
the simulation point of view are generally the rate con- 
stant for the rare transitions between the initial and final 
states and the properties of the Transition Path Ensemble 
(TPE) - the (correctly weighted) collection of transition 
trajectories. When computing these quantities, it is very 
important to know the statistical error in the calculated 
value, and the likely cost of the computation. In this pa- 
per, we derive approximate expressions for these quan- 
tities, for three rare event simulation methods which we 
proposed in a recent publication Q]. These expressions 
turn out to be surprisingly accurate for simulations of a 
model rare event problem. Our results allow us to quan- 
tify the computational efficiency of the three methods. 

The three "FFS-type" simulation methods allow the 
computation of both the rate constant and the transition 
paths for rare events in equilibrium or non-equilibrium 
steady-state systems with stochastic dynamics. In all 
three methods, a series of interfaces are defined between 
the initial and final states. The rate constant is given by 
the flux of trajectories crossing the first interface, mul- 
tiplied by the probability that these trajectories subse- 
quently reach B. The latter probability is computed by 
carrying out a series of "trial" runs between successive in- 
terfaces; this procedure also generates transition paths, 
which are chains of connected successful trial runs. The 
methods differ in the way the trial runs are fired and the 
transition paths are generated. In the "forward flux sam- 
pling" (FFS) method, a collection of points is generated 
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at the first interface and trial runs are used to propagate 
this collection of points to subsequent interfaces - thus 
generating many transition paths simultaneously. In the 
branched growth (BG) method, a single point is gener- 
ated at the first interface and is used as the starting point 
for multiple trial runs to the next interface. Each suc- 
cessful trial generates a starting point for multiple trials 
to the following interface, so that a "branching tree" of 
transition paths is generated. In the Rosenbluth (RB) 
method, a single starting point is chosen at the first in- 
terface, multiple trial runs are carried out, but only one 
successful trial is used to propagate the path to the next 
interface - thus unbranched transition paths are gener- 
ated. In this method, a re-weighting step is needed to 
ensure correctly weighted transition paths. 

A range of simulation techniques for rare events in 
soft condensed matter systems are currently available. 
In Bennett- Chandler- type methods, the rate constant is 
obtained via a computation of a free energy barrier 0. 
In Transition Path Sampling (TPS) §, transition tra- 
jectories (paths) are generated by shooting forwards and 
backwards in time from already existing paths, and are 
then sampled using a Monte Carlo procedure. The rate 
constant is obtained via the computation of a time cor- 
relation function. Bcnnctt-Chandler-type methods and 
TPS are suitable for systems with stochastic or determin- 
istic dynamics, but they require knowledge of the steady 
state phase space density, which means that the system 
must be in equilibrium. While the FFS-type methods are 
only suitable for systems with stochastic dynamics, they 
do not require the phase space density to be known and 
can therefore be used for non-equilibrium steady states 
not satisfying detailed balance. To our knowledge, the 
only other path sampling method that is suitable for non- 
equilibrium systems is that proposed recently by Crooks 
and Chandler 0, which adopts a "TPS" -type methodol- 
ogy, generating new stochastic paths from old paths by 
changing the random number history. 

The origin of the efficiency of the FFS-type methods 
is that they use a series of interfaces in phase space be- 
tween the initial and final states to divide up the tran- 
sition paths into a series of connected "partial paths" . 
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These partial paths are generated in a ratchet-like man- 
ner - i.e. once a particular interface has been reached, 
the system configuration is stored and is used to initi- 
ate trial runs to the next interface. Many other rare 
event techniques also use a series of interfaces in phase 
space. In Transition Interface Sampling (TIS) [5j and 
Partial Path Transition Interface Sampling (PPTIS) @, 
interfaces are used to facilitate the generation of transi- 
tion paths by a TPS-like procedure. In Milestoning [jj, 
trajectories are generated between interfaces assuming a 
steady-state distribution at each interface, while string 
methods 0, use a series of planes in phase space to al- 
low a trajectory connecting the initial and final states to 
relax to the minimum free energy path. The advantages 
of the FFS-type methods over other transition path and 
rate constant calculation methods are that no assump- 
tions are made about "loss of memory" during the transi- 
tion, no a priori knowledge is required of the steady state 
phase space density, and the rate constant is obtained in 
a simple and straightforward way. We have recently be- 
come aware that the BG method bears resemblance to 
the RESTART method, used for simulating telecommu- 
nications networks 0, 0, 113 (this approach was orig- 
inally introduced by Bayes [13| ). The efficiency of that 
method has also been analysed A related method, 
known as Weighted Ensemble Brownian Dynamics, has 
been applied to protein association reactions . 

The key aim of a rare event simulation technique is 
to calculate the rate constant, or in some cases, obtain 
the TPE, with enhanced efficiency compared to brute 
force simulations. However, quantifying the efficiency of 
a particular simulation method is often difficult. Our aim 
in this paper is to derive simple but accurate expressions 
for the computational cost and statistical accuracy of the 
three FFS-type methods. We define the "efficiency" of 
the methods to be the inverse of the product of the cost 
and the variance in the calculated rate constant; our re- 
sults then allow us to analyse the efficiency of the meth- 
ods in a systematic way. From a practical point of view, 
we expect the expressions derived here to be of use to 
those carrying out simulations in two ways. Firstly, when 
faced with a rare event problem, one often has a limited 
amount of computer time available, and specific require- 
ments as to the desired accuracy of the calculated rate 
constant. Analytical expressions for the cost and statis- 
tical accuracy would allow one to estimate, before begin- 
ning the calculation, whether the desired accuracy can be 
obtained within the available time, and thus to make an 
informed decision as to which, if any, method to use. Sec- 
ondly after completing a rate constant calculation, one 
needs to obtain error bars on the resulting value - this is 
especially important for rare events, where both experi- 
mental and simulation results can be highly inaccurate. 
In general, error estimation requires the calculation to 
be repeated several times, which is computationally ex- 
pensive. However, if analytical expressions were available 
for the statistical accuracy in terms of quantities which 
were already measured during the rate constant calcu- 



lation, one could obtain the error bars on the predicted 
rate constant, to within reasonable accuracy, without the 
need for lengthy additional calculations. In this paper, 
we derive such analytical expressions. 

Approximate expressions are derived for the cost, in 
simulation steps, and for the variance in the calculated 
rate constant, for the three FFS-type methods. We ini- 
tially treat the simple case where all trials fired from one 
interface have equal probability of succeeding. We then 
move on to the more realistic case where the probabil- 
ity of reaching the next interface depends on the identity 
of the starting point. To this end, we include in our 
calculations the "landscape variance" - the variance in 
the probability of reaching the next interface, due to the 
characteristic "landscape" for this particular rare event 
problem. Our expressions are functions of user-defined 
parameters, such as the number of trial runs per point at 
a particular interface, as well as parameters characteriz- 
ing the rare event problem itself, such as the probability 
that a trial run succeeds in reaching the next interface. 

We analyse the efficiency of the three methods as a 
function of the parameters, for a "generalized" model 
system. We find that the optimum efficiency is similar 
for all three methods, but that the effects of changing the 
parameter values are very different for the three methods. 
In particular, the BG method performs well only within 
a narrow range of parameter values, while the FFS and 
RB methods are more robust to changes in the param- 
eters. The RB method has consistently lower efficiency, 
due to its requirement for an acceptance/rejection step - 
however, RB may be more suitable for applications where 
analysis of transition paths as well as rates is needed, or 
where storage of configurations is very expensive. 

To test the accuracy of our predictions in the con- 
text of a real simulation problem, we then apply the 
three FFS-type methods to the two-dimensional non- 
equilibrium rare event problem proposed by Maier and 
Stein 0, 0, 0| . We measure the computational cost 
of the methods and the variance in the final value of the 
rate constant, and we compare these to the cost and vari- 
ance predicted by the expressions derived earlier. We 
find that the expressions give remarkably good predic- 
tions, both for the cost and the variance. This suggests 
that the expressions can, indeed, be used to give accurate 
and easy-to-calculate error estimates for real simulation 
problems. 

In Section [HI we briefly describe the three FFS-type 
methods. Expressions for the computational cost and for 
the statistical error in the calculated rate constant are 
derived in Section ITTT1 In Section llVl these expressions 
are shown to be accurate for the two-dimensional non- 
equilibrium rare event problem proposed by Maier and 
Stein. Finally we discuss our conclusions in Section Ivl 
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II. BACKGROUND: FFS-TYPE METHODS 

The FFS-type methods use the "effective positive flux" 
expression for the rate constant, which was rigorously 
derived by Van Erp et al H El H E3 • The rare 
event consists of a transition between two regions of phase 
space A(x) and B(x), where x denotes the coordinates of 
the phase space. The transition occurs much faster than 
the average waiting time in the A state. We assume that 
a parameter X(x) can be defined, such that A < Xa in 
A and A > Ag in B. A series of values of A, A . . . A„, 
are chosen such that Ao = Xa, A„ = Xb and A, < Ai_i. 
These must constitute a series of non-intersecting sur- 
faces in phase space, such that any transition path lead- 
ing from A to B passes through each surface in turn. 
This is illustrated in Figure ^ 




FIG. 1: Schematic illustration of the definition of regions A 
and B and the interfaces Ao . . . X n (Here, n = 3). Three 
transition paths are shown. 



The rate constant Uab can be expressed as [2C 



kAi 



,1.0 



P(X n \X ). 



(1) 



In Eq. pp. Ha is a history-dependent function describ- 
ing whether the system was more recently in A or B: 
hA = 1 if the system was more recently in A than in B, 
and hA — otherwise H flU |20) . The over-bar denotes a 
time average. <&a,j is the flux of trajectories with hA = 1 



that cross A, for the first time 



those trajectories 



that cross Aj, having been in A more recently than any 
previous crossings of Aj. P(Aj|A,-) is the probability that 
a trajectory that comes from A and crosses Xi for the 
first time will subsequently reach Aj before returning to 
A: thus P(X n | Ao) is the probability that a trajectory that 
leaves A and crosses Ao will subsequently reach B before 
returning to A. Eq.Q states that the flux of trajecto- 
ries from A to B can be expressed as the flux leaving A 
and crossing Ao, multiplied by the probability that one 
of these trajectories will subsequently arrive at B rather 
that returning to A. P(A„|Ao) can be expressed as the 
product of the probabilities of reaching each successive 
interface from the previous one, without returning to A: 



P(K\X ) = Y[P(X i+1 \Xi) 



(2) 



i=0 



For simplicity of notation, in what follows, we define 
Pb ^P(A_ Il |A ), Pi = P(A i+ i|Aj), qi ee 1 - pi and 
<3> ee $>A.o/hA- We also use the superscript "e" to in- 
dicate an estimated value of a particular quantity. 

Previously, we described in detail three different ap- 
proaches -the "forward flux sampling" (FFS), "branched 
growth" (BG) and "Rosenbluth" (RB) methods - to cal- 
culating kAB , based on expressions JTJ and J5J jl El ■ 
For completeness, we briefly repeat the description here. 



A. Forward flux sampling 

In FFS, the flux $ is measured using a free simulation 
in the basin of attraction of region A. When the system 
leaves A and crosses Ao for the first time (since leaving 
A), its phase space coordinates are stored and the run is 
continued. In this way, a collection of Nq points at Ao is 
generated, after which the simulation run is terminated. 

The probabilities pi are then estimated using a trial 
run procedure. Beginning with the collection of points at 
Ao, a large number Mq of trials are carried out. For each 
trial, a point is selected at random from the collection 
at Ao. This point is used to initiate a simulation run, 
which is continued until the system either crosses the 
next interface Ai, or re-enters A. If Ai is reached, the 
final point of the run is stored in a new collection. After 

Mo trials, po is given by 7Vi°' /Mo, where Ai ' is the 
number of trials which reached X±. The probability p\ 
is then estimated in the same way: the new collection of 
points at Ai is used to initiate M\ trial runs to A2 (or 
back to A), generating a new collection of points at A2, 
and so on. Finally, the rate constant is obtained using 
Eqs CQ) and ©. 

FFS generates transition paths according to their cor- 
rect weights in the TPE [j],|2l|. In order to analyse these 
transition paths, one begins with the collection of trial 
runs which arrive at Xb from A„~i and traces back the 
sequence of connected partial paths which link them to 
region A. The resulting transition paths are branched 
- i.e. a single point at Ao can be the starting point of 
multiple transition paths. 



B. The branched growth method 

In the BG method, which was inspired by techniques 
for polymer sampling 0, [2^, [2^, branched transition 
paths are generated one by one, rather than simultane- 
ously, as in FFS. The generation of each path begins with 
a single point at Ao, obtained using a simulation in the 
basin of attraction of A, as in the FFS method. This 
point is used to initiate fco trial runs, which are contin- 
ued until they either reach Ai, or return to A. Each of 

the Ng°^ end points at Ai becomes a starting point for ki 
trial runs to A2 or back to A. Each of the Ng 1 ^ successful 
trial runs to A2 initiates ki trials to A3, and so on until A„ 
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is reached. An estimate Pg of Pb is obtained as the total 
number of branches that eventually reach A„ , divided by 
the total possible number: P B — ivj™ -1 "'/ ]Xi=o ^i- If, at 
any interface, no trials were successful, P% = 0. To gen- 
erate the next branching path, we obtain a new starting 
point at Ao from the simulation in the basin of attraction 
of A. After many branching paths have been generated, 
an average is taken over the PJj values of all the paths. 
The flux $ is meanwhile obtained from the simulation 
run in region A. The branched transition paths that 
are generated in the BG method are correctly weighted 
members of the TPE We note that the BG method 
bears resemblance to methods developed for telecommu- 
nication networks flfj [ [ill and to a method used for 
protein association [Tj|. 

C. The Rosenbluth method 

The RB path sampling method is related to the 
Rosenbluth scheme for sampling polymer configurations 
0, l24l l25| . The RB method generates unbranched tran- 
sition paths, one at a time. An initial point at A is 
obtained using a simulation in the A basin, which is con- 
tinued until the trajectory crosses Ao for the first time, 
as in the FFS and BG methods. This point is used to 
initiate ko trials, which are continued until they either 
reach Ai or return to A. If > of these trials reach 
Ai, one successful trial is selected at random and its end 
point at Ai is used to initiate k\ trials to A2 or back to 
A. Once again a successful trial is chosen at random and 
the process is repeated until either no trials are success- 
ful or A„ is reached. The generation of the next path 
then begins with a new point at Ao, obtained using the 
simulation run in the A basin. 

The Rosenbluth method as outlined above does not, 
however, generate paths according to their correct 
weights in the TPE: for correct sampling, paths must be 
re- weighted by a "Rosenbluth factor". The Rosenbluth 
factor for a partial path up to interface i is given by: 

i-l 

Wi = l[NP (3) 

Note that the re-weighting factor Wi depends on the 
number of successful trials obtained at all the previous 
interfaces, while generating the path up to Ai. The cor- 
rect re- weighting can be achieved using a Metropolis-type 
acceptance/rejection scheme 0, in which a newly gener- 
ated path is either accepted or rejected based on a com- 
parison of its Rosenbluth factor with that of a previously 
generated path. Ensemble averages of any quantity of 
interest are then taken over all accepted paths. Here, the 
quantity which we wish to calculate is the probability pi 
that a trial run fired from Ai will reach Aj+i, for each 
interface i. When we fire ki trial runs from Ai, we obtain 
an estimate for pc pi = Ns /ki. We require the cor- 
rectly weighted ensemble average for pi at each interface 



i; we note, however, that the same procedure could also 
be used to calculate the ensemble average of any other 
property of the ensemble of paths from Ao to Ai. 

From a practical point of view, each interface has as- 
sociated with it two values of Wi and p\. The first set 

of values: and p\ , are associated with the tran- 

sition path that is currently being generated (the "new" 

path). depends on the number of successful trials 

generated in creating this transition path as far as Aj, 
and p^ = Ns /ki depends on the number of successful 
trials fired from the point at Ai to Aj+i. The other set of 

values, andp* , are the "old" values for this inter- 
face. These values correspond to the last "acceptance" 
event at this interface. 

The recipe for obtaining k^B within the RB method is 
as follows. Transition paths are generated as described 
above. When the path generation procedure reaches Aj, 
we calculate the Rosenbluth factor (using Eq.@) 

and we fire ki trial runs to obtain p| = Ng /ki. We 
then calculate the ratio /W^ and draw a random 

number < s < 1. If s < W^/wj°\ an acceptance 
event takes place. In this case, the previous values of 
and pi are replaced by the newly obtained values 
W<; n) and pf n) . If, however, s > wj; n) /W$ o) , a rejec- 
tion occurs and W$ and p^ remain unchanged for 
this interface. Regardless of the outcome of the accep- 
tance/rejection step, the accumulator for the probability 

p\ is incremented by the current value of pf - this may 
be either a newly generated value (if an acceptance just 
occurred) or an old value that may have been already 
added to the accumulator several times (if several rejec- 
tions have happened in a row). To proceed to the next 
interface, a successful trial run is chosen out of those that 
have been newly generated, and its end point at Aj+i is 
used as the starting point for ki+\ trial runs to Aj+2- 
A corresponding acceptance/rejection step is then car- 
ried out at Aj+i. We note that the "old" values W^ 

and j?i^ for different interfaces need not correspond to 
the same transition path. After many complete tran- 
sition paths have been generated, kAB is obtained using 
Eq. (Q , where an estimate of the flux $ is calculated from 
the simulation run in region A. A "pseudo-code" corre- 
sponding to the above procedure is given in our previous 
publication [l|, together with a description of an alter- 
native, "Waste Recycling" [2(| re-weighting scheme. In 
this paper, however, we shall consider only the Metropo- 
lis acceptance/rejection approach. 



III. COMPUTATIONAL EFFICIENCY 

In this section, we derive approximate expressions for 
the computational efficiency of the three methods. Fol- 
lowing Mooij and Frenkel |27|. we use the following defi- 
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nition for the efficiency, £ : 



£ = cv 



(4) 



In Eq. ©, C represents the computational cost, which 
we define to be the average number of simulation steps, 
per initial point at Ao- The statistical error in the esti- 
mated value k AB of the rate constant is represented by 
V. Denoting the mean (expectation value) of variable u 
by E[u] and its variance by V[u], we define V to be the 
variance V[fc^ B ], per initial point at Ao, divided by the 
square of the expectation value E[k e AB \: 



V = 



N V[k 



AB\ 



(E[k 



V\k 



AB\ 



AB\ 



k 2 
K AB 



(5) 



where iVo is the number of starting points at Ao used in 
obtaining the estimate k AB . The expectation value of 
k e AB is, of course, the true rate constant: E[k e AB ] = k AB . 
The error bar for k AB is given by kAB y/V/No. 



A. Computational Cost 

We define the computational cost C of a particular 
method to be the average number of simulation steps 
required by that method, per starting point at Ao- In 
making this definition, we ignore any other contributions 
to the CPU time, such as memory storage. To estimate 
the value of C, we consider a generic system that makes 
a rare transition between states A and B. A parameter 
A and interfaces Ao . . . A„ are chosen as in Section ITT1 

There are two contributions to the cost C. The first 
is the average cost R, in simulation steps, of generating 
one starting point at Ao- This is related to the flux $ 
from the A region to Ao by R = l/($c£i), where dt is the 
simulation timestep. 

The second contribution to C is the cost of the trial 
run procedure. We first consider the cost Ci of firing one 
trial run from interface Aj. The run is continued until it 
reaches either the next interface Aj+i (with probability 
Pi), or the boundary Xa of region A (with probability 
qi). We make the assumption that the average length 
(in simulation steps) of a trajectory from interface Aj to 
another interface Xj is linearly proportional to |Aj — Aj|, 
with proportionality constant S. Ci is then given by: 



C 4 = S [pi(A i+ i - Xi) + q t {Xi - X A )\ 



(6) 



The basis for the assumption of linearity in Eq. © is that 
we suppose that the system undergoes one-dimensional 
diffusion along the A coordinate in the presence of a "drift 
force" of fixed magnitude. For an equilibrium system, 
the origin of the drift force is the free energy barrier. 
Farkas and Fiilop have presented analytical solutions [2^ 
for the mean time to capture for a particle undergoing 
one-dimensional diffusion with constant drift force, in the 
presence of two absorbing boundaries. In Appendix ^ 



we show how these results lead to Eq.©. Eq.© is shown 
to be valid for the two-dimensional Maier-Stein problem 
in Section ITVAl ( Figure I7|l . 



Expressions for the cost 

Given Eq.©, we can compute the average cost C per 
starting point at Ao of the three methods. 

In FFS, we make Mj trial runs from interface i and, 
providing at least one of these is successful, we proceed 
to the next interface i + 1. In practice, M, is expected to 
be large enough that at least one trial run reaches Ai+j.. 
In this case, the expected cost per starting point at Aq 
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C m = R+ — Y j M i C i 



(7) 



Defining ki such that ki = Mi/No, Eq.Ql can be rewrit- 
ten as: 



-FffB 



R 



n-l 

E 

i=0 



(8) 



If, however, Mj is small, we must take account of the 
possibility that none of the trial runs from A; reach Aj+i. 
In this case, the FFS procedure is terminated at interface 
i and the cost is accordingly reduced. Since the proba- 
bility of reaching interface i > is Y[]Jo fl ~ ^j^) (this 
is the probability that at least one trial is successful at 
all interfaces j < i), Eq.© is replaced by: 



-Fffs 



R + fcoCo 



n-l 

•E 

i=l 



i-1 

kiCi 

3=0 



1-q 



Nok, 



(9) 



Although the cost is reduced by failing to reach later in- 
terfaces, this of course results in a less accurate prediction 
of the rate constant, since the terminated FFS calcula- 
tion makes no contribution to the estimate of pi for later 
interfaces. This will be reflected in our expression for the 
statistical error in Section UlI Bl 

We now turn to the BG method. Here, we generate 
a "branching tree" of paths, with N% points at interface 
i originating from a single point at Ao- We fire ki trial 
runs for each of these Ni points. The average value of Ni 
is: 



(i>0) 



(10) 



Of course Nq — 1. The average cost per starting point 
at Aq is therefore: 



n-l 



C bg = R+Y^hCiNi 



(11) 



i=0 



R + k C + Y 



kiCi J J /',/•', 

3=0 







Finally, we come to the RB method. In this algorithm, 
we generate unbranched paths by firing ki trials from in- 
terface i, choosing one successful trial at random and pro- 
ceeding to interface i + 1. If no trial runs are successful, 
we start again with a new point at Ao- The probability 

of reaching interface i > is ]TV 1 (l-q^y The cost 

of the RB method, per starting point at Ao, is therefore: 



a cost that increases dramatically with increasing k or 
n. This effect is due to the rapidly increasing number of 
branches per starting point at Ao- In this regime, the FFS 
and RB methods converge to the same cost, since Eqs © 
and (|12fl become equivalent when 1 — q k w 1 — q N ° k w 1. 

B. Statistical Error 



C rb = R + k C + 



3=0 



(12) 



Once again, the "price" of failing to reach later interfaces 
will be paid in the form of an increased variance in the 
calculated rate constant. The effect of the Metropolis 
acceptance/rejection step in the RB method appears only 
in the variance in k e AB fSection llll Bp . and not in the cost. 



Illustration 
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FIG. 2: Cost C/R, for evenly spaced interfaces, pi — p, 
ki = k, R = S, No = 1000 and P B = 10~ 8 . (a): C/R as a 
function of k, for n = 5. (b): C/R as a function of n, for 
k = 25. 



For the purposes of illustration, let us consider a hypo- 
thetical rare event problem for which Ao = A^ = and 
A n = Ab = 1< We suppose that the interfaces are evenly 
spaced in A, have equal values of pi, and that the firing 
parameter ki is the same at each interface: i.e. Xi = i/n, 

Pi = P B n (from Eq.Q) and ki — k. We also suppose 
that R = S and Nq = 1000. The resulting values of the 
cost C, obtained from Eqs JJjJ, l|ll|) and (|12f> . are plotted 
in Figure^ and b as functions of k and n. In the regime 
of small k or small n (implying small p) , the BG and RB 
methods converge, while the cost of the FFS method is 
higher. This is because, for BG and RB, the probability 
of reaching later interfaces is low and the cost is domi- 
nated by the trial runs fired from early interfaces. The 
FFS procedure is less likely to be terminated at early in- 



terfaces (note the factor of 1 — qi -i in Eq. @ as opposed 
to 1 — q k ' in Eq. (|12[0 . and is therefore more expensive, 
per initial point at Ao . In the regime of large k or large 
n (implying large p), a different scenario emerges. Here, 
the BG method becomes by far the most expensive, with 



We now turn to the relative variance V in the estimated 
value k e AB of the rate constant, per starting point at Ao- 
k AB is the product of the estimated flux through Ao, 
multiplied by the estimated probability of subsequently 
reaching B: k e AB = <5> e P B (Eq.Q). 

In this paper, we shall ignore the error in $ e . <I> e is 
obtained by carrying out a simulation run in the basin 
of attraction of A and measuring the average number of 
simulation steps between successive crossings of Ao (com- 
ing directly from A). As long as Ao is positioned close 
enough to the A region, the simulation run in A can 
be made long enough to estimate $ with high accuracy, 
with a computational cost that is minimal compared to 
the cost of estimating P B . We therefore obtain: 



V = N Q - 



V[k 



ABl 



N 



3> 2 V[P B ] 



N 



r B 



(13) 



{E[k AB \y "{mp B \) 2 

In Ea. (|13|) . we have used the general relation [2^| 

V[ax] = a 2 V[x] (14) 

where a is a constant. 

In what follows, we shall make the important assump- 
tion that the numbers N^ of successful trial runs at dif- 
ferent interfaces i are uncorrelated - i.e. that if, during 
the generation of a transition path, one is particularly 
successful or unsuccessful at interface i, this will have no 
effect on the chances of success at interface i + 1. In 
reality, of course, there will be correlation between inter- 
faces, especially if the interfaces are closely spaced or the 
system dynamics have a large degree of "memory" . We 
expect this assumption to be the major limiting factor in 
the applicability of our results to real systems; however, 
as we shall see in Section Hvl the results are surprisingly 
accurate for the two-dimensional Maier-Stein problem. 
We expect that the expressions derived here could be 
modified to include the effects of correlations between 
interfaces; for highly correlated systems this may prove 
necessary. 

Expressions for the variance 

The basis of our analysis is the fact that on firing ki 
trial runs from interface i, the number of successful trials 



Ns is binomially distributed |29|, with mean 
E[N^} = hpi 



(15) 
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and variance 



V[N&] = k iPiqi 



(16) 



For now, we assume that all trial runs fired from interface 
Xi have equal probability p t of reaching Aj+i. This as- 
sumption will later be relaxed. We shall need to express 
the variance in P B in terms of the variance V\pf] in the 
estimated values of Pi at each interface. To do this, we 



recall that P% = Ui=qJ. 
the following relation [2£ 



(Eq.@), and we make use of 



V[f(x,y, ...)] = 



dx 



V[a 



df 
dy 



V[y] 



(17) 



where f(x, y, . . . ) is a function of multiple uncorrelated 
variables x,y, . . . and the partial derivatives are evalu- 
ated with all variables at their mean values. By "uncorre- 
lated variables" we mean that the covariance Cov[m, v] = 
for all pairs of variables u and v. Identifying x, y . . . 
with pf,pf +1 ... and taking f(jp e . . .p e n _ x ) = T\7=o Pi, we 
find that df/dpf = [n"=o Pjl/Pi = PbM* 80 that 



V\Pi 



J2 E 

i=0 



pi 
Pi 



v[p1 



p2 



E 

i=l 



V[p! 



Pi 



(18) 



We now use the above results to calculate V for the 
FFS method. In this method, we begin with a collection 
of N points at A . For each interface, p\ is obtained 
by firing Mj = Noh trial runs: pf = Ng /Mi, where 
is the number of trials which reach Ai+i. Using 
V[N^ ] ]/Mf. Using Eq.lJTSJl, we find 
MiPiqi. Noting also that E[pf\ = pi and 



Eq.lHU, V\p\ 

that V[N< 
using Eq. lfTBl) . we obtain 



Wi 



U P * M * N ° <=o Plk 



JL (19) 



and from Eq.JTSJ 



V 



ffs 



n-l 

E 

i=0 



Piki 



(20) 



As for the cost calculation, we have assumed that Mj 
is large enough that there is always at least one trial 
run which reaches the next interface. If this is not the 
case, we must also take account of the possibility that 
interfaces i > may not be reached. The probability of 

reaching interface i > is Il}=o fl — ftT^)' so 



( i - * 



Mi 



1 



% ) 



(21) 



Eq. (|21|l is written in this form so that for i = 0, we 
recover V[p{|] = PiQi/Mi. Eqs (fTT?f) and (|2*n|l must then 



be replaced by: 



E 



(22) 



and 



n-l 



V 



ri's 



E ~ 



i=0 



1-9. 



ifi=o ( 



N„k, 



(23) 



We now turn to the BG method. Here, we begin with 
a single point at Ao- From this point, we generate a 
branching "tree" of paths connecting A to B. The value 
of Pb is estimated by 



P' 



N, 



(n-l) 



Yii=0 hi 



(24) 



where ivi n ^ is the total number of trials reaching 
A„ = Xb- We denote the number of points in the 
branching tree at interface i by iVj. For a given num- 
ber N n _i of points at A n _i, the total number of tri- 
als fired is N n —ik n -i and the variance in Ng ^ is 

V r [JV, (n ~ 1) |JV fl _i] = JV„-ifc»-ii>»-ig„-i (using Eq.©). 
However, the situation is complicated by the fact that 
N n -i itself varies; in fact, N n _i is simply the number 
of successful trial runs reaching A„_i from A ra _a, and in 
general: 



N 4 = 



>0] 



(25) 



At this point, we need to calculate the variance in a 
quantity Y which is conditional upon the value of another 
quantity X. Here, and several times in the rest of the 
paper, we will use the general relation 



V[Y] = E [V[Y\X]} + V [E[Y\X]} 



(26) 



where the mean and variance on the r.h.s. of Ea. H26|) 
are taken over the distribution of values of X. Since 

B[JVi B - 1) |JV n _ 1 ] = Nn-ihrtr-iPn-l, 



V 



E 



N^\N n ^ = k 2 n _-j? n _^V [JV n _i] (27) 



kn-lPn-lV 



(using Eqs ifTljl and We also know that 



E 



V 



n-2 

= k n -ip n -iq n -i Y[ k iPi 



so that 



n-l 



v{Ni n -^} = n k iPi + ki_ lP i^v 

j=0 



i=0 



(29) 
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Using the same arguments, we can generalize Ea. H29(l to 

^n^^+^M^ 1 '] [i>o](so) 

q i k l p l [i = 0] 

Using Eq. i|30[l , we can solve Eq. I|29|l recursively, to obtain 
V[Ns n ~^}. Using Eqs.I23} and @2J, we then arrive at 
the variance in the estimated value of Pb'- 



(31) 



where we have divided by No to account for the fact that 
P B is calculated by averaging results over N starting 
points at Ao- We then obtain from Ea. i|13|) : 



v bs = E 



i=0 



(32) 



Finally, let us derive the equivalent expression for the 
RB method. Here, we again use Eq. (|18() . If we ignore for 
the moment the effect of the acceptance rejection step, 
we can use Eas. (|16|) and (|14|l to obtain an expression for 
the variance in pf : 



v Th \pf 



Vili 



Noki ii; ,(i % 



(33) 



where we have taken account of the fact that the proba- 
bility of reaching interface i > is n}=o (-*- — % y > an< ^ 
that the pf value is averaged over N separate path gener- 
ations. Eq.(j23 is very similar to the FFS result, Ea.l|2"T]l. 

The Metropolis acceptance/rejection step (described 
in Section ^) increases the variance in pf. On reach- 
ing interface i, we fire ki trials and obtain an estimate 
^e,(n) _ jyiv e j|; ner acce pt or reject this esti- 
mate. If we reject, makes no contribution to the 
average value of pf - instead, the previously accepted 
estimate, J>\'^°\ is added to the average, even though 
pl'° was already added to the average in the previous ac- 
ceptance/rejection step. If, instead, we accept , it 
makes a contribution to pf, and, if the subsequent esti- 
mates happen to be rejected, it may repeat this contribu- 
tion multiple times. The final estimate, pf, is therefore 
an average over all the values of N^ /ki that were gen- 
erated, weighted by the number of times Q that each of 
these values contributed to pf: 



N^/h 



N, 



(34) 



where the sum is over all generated Ng /ki values and 



Ng is the total number of these. In fact, 



N® = N 



n;= ( 



=n I"? 



1_n 
-Qt 



(35) 



since the number of times we fire trials from A,; is simply 
the number of times we begin a path generation from Ao 
and succeed in reaching A^. Using Eq.l|14|). the variance 
pf is then 



2 V 



N^/h 



V 



N 



(01 



Of 



(36) 

(assuming that the distributions of the stochastic vari- 
ables Qi and [Ns^ji are uncorrelated) . Ea. (|36|l is equiv- 
alent to: 



V 



V[ P f 



N, 



k?NP 



E Q 2p (® 



PiQi 



2=o 



k,N, 



(i) 



E Q 2p (Q) 



2=0 



(37) 

In order to find the distribution P(Q), we define a new 
variable 6>; is the probability that we accept a newly 
generated estimate = Ns /ki. P{Q) is then: 

P{Q) 
p{Q) 



(i 



Q = 
Q > 



(38) 



Ea. (l38|) can be understood as follows: Q — corresponds 
to a pl'^ value that is generated but is immediately re- 
jected and therefore contributes zero times to the aver- 
age. This occurs with probability 1 — 0j. Q > corre- 
sponds to a pf^ value that is generated and accepted 
(with probability 9i) - the next Q— 1 values that are gen- 
erated are rejected (with probability (1 — l^)^ -1 ), then 
finally a new value is generated which is accepted (with 
probability Oj), so that the original value ceases to con- 
tribute to the average. The distribution has the 
property that |3(j 



E Q 2p (Q) = 

Q=0 



2-i 



(39) 



so that Eq. (|37|l for the variance in pf per point at Xi 
becomes 



V[ P f] = 



PiQi 



kN 

Using Ea. H35|) . we obtain: 

"(2-. 



(i) 



(40) 



V Tb \pf] = 



Pi<li 
N ki 



1 n 
1 ~1i 



lfi=o(l-#) 



(41) 



Comparing to Ea. (|33|) . we see that the effect of the ac- 
ceptance/rejection step is to multiply V\pf\ by a factor 
(2 — 6i)/9i. Using Ea. (|18|) . the relative variance in P B is: 



V^[P B ] 



ft (2 



p2 
B 



N pih 



n;= ( 



=o 1-9/ 



(42) 



so that using Ea. lfT^l . 

q l (2-0*) 



V 



rb 



E 



1-3? 



Piki 



(43) 



We show in Appendix [E] that the acceptance probability 
#i for i > [note that 9q = 1] can be approximated as: 



1 _VtF 

2 4 



2E,-f ( -± 
2 



(i > 0) (44) 

where Erf (a;) is the error function: Erf (a;) = 
(2/Vtt) Jo e't dt, and (Tj is given by: 



3=0 



(l-g 3 ')fr- 

kjpj 



(45) 



Eqs and can be substituted into Eq. (|4^|) to give 
a complete expression for the relative variance in the es- 
timated rate constant for the RB method. 



Illustration 
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FIG. 3: Relative variance V, for pi = p, ki = k and Pb = 
10 -8 . The circles show the function Y17=o Qi/iPi^i)- ( a ) : ^ 
as a function of k, for n = 5. (b): V as a function of n, for 
k = 25. 



Returning to the hypothetical rare event problem 
with evenly spaced interfaces introduced above, Figure 
shows V function of k (for n = 5) and of n 



(for fc = 25), for pi — p = P 



l/ti 



fc, 7V 



1000 



and Pb = I0~ 8 . The circles show the limiting form 
Y^i=o Qi/(Piki), which is in good agreement with the FFS 
results, since 1 — q N ° k « 1. For small k or small n (small 
p), the RB and BG results tend to converge, since the 
probability of reaching later interfaces is small and the 
results are dominated by the early interfaces. In this 
regime, the FFS method gives the smallest variance, since 
the chance of terminating the trial run procedure at early 
interfaces is lower than for the other methods. 

It is interesting to compare expressions (12311 . and 
(|43[l . All three expressions are of the form 



v = E 



£=0 



(46) 



However, Xi takes different forms for the three methods: 



X 



lis 



j=o0--9j ) 



(1 



.7=0 



and 



X" h = 



Bi nl=o(i 



qf 



(2-0;) 



(47) 



(48) 



(49) 



We note that Xf s > X 4 rb , so that V ffs is always less 
than V rb , even for 4 = 1. Both Xf s and A[ b are al- 
ways less than unity: V approaches the limiting form 
SiLo Qi/(piki) f rom above as fej increases (in fact in 
Fig. it takes this form for all k) and V rb approaches 
Si^To — Oi)qi/ {jpikiOi). For the BG method, however, 
A bg can increase indefinitely as fcj increases, so that this 
method produces the smallest variance for large fcj, as 
in Figure OK. However, comparing with Figure we see 
that this is also the regime in which the BG method be- 
comes very expensive. 



Landscape Variance 

So far in our analysis, we have assumed that all the 
points at interface Xi have to same pi value - i. e. that on 
firing a trial run to A-j+i we have the same probability of 
success, no matter which point at A^ we start from. In 
reality, this is not the case; we expect there to be a distri- 
bution of pi values among the points at each interface Xi. 
We call the variance of this distribution the "landscape 
variance" Ui at interface i, and we expect it to make a 
contribution to the variance in P|j. We now extend our 
analysis to include the potentially important effect of the 
landscape variance. 

Let us suppose that each point j at A^ has an associated 
probability p[ that a trial run fired from that point will 
reach Aj+i. The distribution of p^p values encountered 
during the rate constant calculation has mean i?[p^] = 

Pi and variance V\p\ ] = Ui. Of course, the values of Ui 
depend on the number and placement of the interfaces. 

In Appendix O we re-derive expressions for the rela- 
tive variance in the estimated rate constant, taking into 
account the landscape variance. The final results are: 



V 
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E 

i=0 
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N ki 



1i 



N k, 
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=o 1 - 9 



Nokj 
J 



(50) 



1U 



where JVj = A/ofej-iPi-i for i > and JV^ = iVo for i = 0. 



ra-1 



v bs = E 



i=0 



and 



V r 



E 

i=0 



(51) 



(52) 



(2 



1 ™ 



the relative variance V (as in Figure calculated with 
U = hp 2 /n (upper curves), U = p 2 /n (middle curves) and 
U = (lower curves). Although the landscape variance 
does not change the general trend that V decreases as k or 
N increases, it does have the qualitative effect that V no 
longer tends to zero (as discussed above). Depending on 
the value of U, the quantitative effects of the landscape 
contribution can be very significant, especially as k or N 
becomes large. 



C. Efficiency 



Comparing Eqs l|5U|) . I|51|) and l|52[) to their equivalent 
forms without landscape variance, l|23|) . I|32l) and (|43H . 

we see that for each interface the "binomial" terms of the 
form piqi/ki are now supplemented by additional terms 
describing the landscape variance. In the limit of very 
large ki, the relative variance no longer tends to zero. In- 
stead, as ki — ► oo (for all i), the FFS and BG expressions 
(|50|l and l|51l) tend to the constant value Uq/pq, while 

the RB expression (E2> tends to J2t=o U i/Pi- while the 
"binomial" contribution to the variance can be reduced 
by firing many trial runs per point, the "landscape" con- 
tribution can only be reduced by sampling many points. 
In the FFS and BG methods, branching paths are gener- 
ated. For very large ki, each point at Ao generates many 
points at subsequent interfaces, so that only Uq remains 
in Eqs 150JI and (|51|l as ki — > oo. In the RB method, 
however, paths are not branched, so that each point at 
Ao corresponds to one (or less than one) point at each 
subsequent interface. In this case, as ki — > oo, all the Ui 
values continue to contribute to V. 
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FIG. 4: Relative variance V in k e AB , as predicted by Eqs 1501 . 
151H and 152H , for the model problem of Figs |5] and |3] with 
Pb = 10 -8 and Ui — U. The upper curves in each group 
correspond to U = 5p 2 /n, the middle curves to U = p 2 /n and 
the lower curves to U = 0. (a): V as a function of k, keeping 
n = 5. (b): V as a function of n, keeping k = 25. 



In Figure 21 w e revisit the simple model problem of 
Figs [3 and |3 adding in the effects of landscape variance. 
We take Ui to be the same for all interfaces: Ui — U . We 
choose, somewhat arbitrarily, U — p 2 /n or U = hp 2 jn. 
These turn out to be quite realistic values for the Maier- 
Stein system discussed in Section IIVI Figure 0] shows 
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FIG. 5: Efficiency £, calculated using Eq.Q, for the simple 
model of Figs|5|01and2] For each method, results are plotted 
with U = p 2 /n (lower curves) and (7 = (upper curves), (a): 
£ as a function of k for n — 5. (b): £ as a function of n for 
k = 25. 



Having calculated the computational cost and the sta- 
tistical accuracy of the three methods, we are now in a 
position to assess their overall computational efficiency, 
as defined by Eq.QJ. Figure |3] shows the efficiency of the 
three methods as a function of k (Fig. |5Ji) and of n (Fig. 
|SJj), for the simple model case of Figs |21 El and 0] Note 
the altered scale on the n axis in comparison to Figures 
121 and For each method, the upper curve shows the 
results without the landscape contribution to the vari- 
ance (U = 0) and the lower curve includes a landscape 
contribution of U — p 2 /n. 

Firstly, we note that the optimum values of £ are of 
the same order of magnitude for all three methods, al- 
though £ is consistently lower for RB, due to the ac- 
ceptance/rejection step. However, the dependence of the 
efficiency on the parameter values k and n is very differ- 
ent for the three methods. For the BG method, the effi- 
ciency shows a pronounced peak, both as a function of k 
and of n. Although for an optimum choice of parameters, 
this method can be the most efficient, its performance is 
highly sensitive to the choice of parameters, decreasing 
sharply for non-optimal values of k or n. The FFS and 
RB methods are much less parameter-sensitive - in fact, 
as long as k or n is not too small, the choice of parame- 
ters appears not to be at all critical for these methods. In 
general, FiglSJseems to indicate that the method of choice 
is FFS, since this method is highly robust to changes in 
the parameters, is the most efficient method at small k or 
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n, and remains efficient as k and n become large. How- 
ever, this interpretation must be treated with care, since 
several important factors are not included in the analysis 
leading to FigEI Firstly, our analysis does not include 
the effects of correlations between interfaces. This has 
the effect that neither the FFS or RB methods shows a 
maximum in efficiency as a function of n in Fig Eh- In our 
simple model, one can always gain more information by 
sampling at more closely spaced interfaces - however, in 
reality, correlations between interfaces are likely to make 
very closely spaced interfaces computationally inefficient. 
Another important factor to be considered is the fact 
that both the FFS and BG methods generate branched 
transition paths. In FFS, in fact, an effect analogous 
to "genetic drift" means that if the number of points in 
the collections at the interfaces is small enough to be of 
the order of the number of interfaces, then all the paths 
that finally reach B can be expected to originate from a 
small number of initial points at Ao If there is "memory 
loss" - i.e. no correlations between interfaces, this may 
be unimportant. However, if the history of the paths 
is important, then the RB method may be the method 
of choice, since this generates independent, unbranched 
paths. Furthermore, the RB method requires much less 
storage of system configurations than FFS (for which a 
whole collection of points must be stored in memory at 
each interface) - for some systems, this may be a signifi- 
cant factor in the computational cost. 

Figure [5] also shows the effects of landscape variance 
on the efficiency of the three methods. Including land- 
scape variance always decreases the efficiency, but pro- 
duces rather few qualitative effects for this simple model 
problem. It is interesting to note, however, that in Figure 
EJi both the FFS and RB methods show a maximum in 
efficiency as a function of k only when the landscape con- 
tribution is included. When the landscape contribution 
is not considered, the equations predict that arbitrarily 
high accuracy can be obtained by firing an infinitely large 
number of trials from a single point. In this example, we 
took the landscape variance to be the same for all in- 
terfaces: Ui — U . However, one can easily imagine that 
for some systems, there is much greater variation among 
transition paths when they are close to the A basin, while 
for others, paths tend to diverge as they approach B. In 
the former case, we can expect the RB and BG methods 
to have an advantage relative to FFS, because in these 
methods, relatively more points are sampled at early in- 
terfaces (since the probability of failing to complete a 
transition path is higher). Conversely, if the landscape 
variance is very large close to the B basin, the BG method 
may be advantageous, since it samples many points at 
later interfaces due to its branching tree of paths. 



our test case, we simulate the two-dimensional non- 
equilibrium rare event problem proposed by Maier and 
Stein 0,0,0. This system has been extensively stud- 
ied both theoretically and experimentally 0, 0, fl^.l3ll 
|3^ | and was also used by Crooks and Chandler Q as 
a test case for their non-equilibrium rare event method. 
We hope that the conclusions obtained for this system 
will also prove to be applicable to more computationally 
intensive rare event problems. 
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FIG. 6: Typical trajectory for a brute-force simulation of the 
Maier-Stein system, with a = 6.67, /i = 2 and e — 0.1. 

The Maier-Stein system consists of a single particle 
moving with over-damped Langevin dynamics in a two- 
dimensional force field. The position vector (xi,X2) of 
the particle satisfies the stochastic differential equation: 

Xi = fi(x)+ti(t) (53) 

where the force field f = (/i, f%) is given by: 

f = (xi — x\ — ax\x\ , —11x2(1 + (54) 

and the stochastic force £ = (^1,^2) satisfies: 

&(*)) = ! m + r)&(i)> =eS(t- r)5 %3 (55) 

This system is bistable, with stable points at (±1, 0) and 
a saddle point at (0,0). If a ^ /j, the force field f can- 
not be expressed as the gradient of a potential. In this 
case, the system is intrinsically out of equilibrium and 
does not satisfy detailed balance. The parameter e con- 
trols the magnitude of the stochastic force acting on the 
particle. For e > 0, the system makes stochastic tran- 
sitions between the two stable states, at a rate which 
decreases as e decreases. Figure El shows a typical trajec- 
tory generated by a brute-force simulation. Here, and in 
the rest of this Section, we use a = 6.67, // = 2.0 (fol- 
lowing Crooks and Chandler Q) and e = 0.1. Eg. 153(1 is 
integrated numerically with timestep 8t = 0.02 [33|. For 
our calculations using the FFS- type methods, we define 
A(x) = x\, \a = Aq = —0.7 and As = A„ = 0.7. 



IV. THE MAIER-STEIN SYSTEM 



A. Measuring the parameters 



In this section, we test the expressions derived in Sec- 
tion IIIII for a real rare event simulation problem. As 



In order to test the expressions of Section IITTI we must 
measure the cost parameters R and S, the probability 
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Pb of reaching B and, for a given set of n interfaces, 
the probabilities {pi} and the landscape variance values 
{Ui}. For most of our calculations, we used n = 7, and 
the interfaces were positioned as listed in Table |U For 
the results of Figs 0>, and lllb . where n was varied, 
we kept the interfaces evenly spaced between Ao = —0.7 
and A n = 0.7. R, the cost of generating an initial point 
at Ao, was measured using a simulation in region A to 
be R — 590 ± 50 steps. In these calculations, points at 
Ao were collected upon every 10th crossing of Ao from A. 
To measure S (the proportionality constant in Eq.@|), 
we carried out an FFS run, measuring the average length 
(in simulation steps) of successful and unsuccessful trials 
from each interface. The results are shown in Figure 
Here, the filled circles show the average length, in simu- 
lation steps, of successful trials from interface Xi (plotted 
on the x axis) to A^+i = Xi + 0.2. Since |Ai — Aj| = 0.2 
for all these trials, Eq.JSJ) predicts that all the filled cir- 
cles should have show the same average trial length. The 
open circles show the average length of unsuccessful tri- 
als, which begin at Ai and end at Xa = —0.7, so that 
|A; — Xj\ = Xi + 0.7: Eq.© predicts that all the open 
circles should lie on a straight line. Combining all the 
data, we obtain an average value of S = 131 steps. This 
value is used to plot the solid lines in Figure The very 
good agreement that is observed between the solid lines 
and the circles implies that the drift-diffusion approxi- 
mation, Eq.©, is reasonable for this problem. The most 
significant deviation occurs for the successful trial runs 
between A = —0.7 and A = — 0.5; these are unexpectedly 
short, perhaps because the "drift force" is weaker in this 
region. 



Interface 

1 
2 
3 
4 
5 
6 



Xi 
-0.7 
-0.5 
-0.3 
-0.1 
0.1 
0.3 
0.5 



Pi 

0.1144 ± 0.0001 
0.2651 ± 0.0002 
0.3834 ± 0.0002 
0.5633 ± 0.0003 
0.7702 ± 0.0003 
0.9152 ± 0.0002 
0.9747 ± 0.0002 



Ui 

0.00350 ± 0.00003 
0.00368 ± 0.00008 
0.0031 ± 0.0003 
0.0021 ± 0.0002 
0.0008 ± 0.0001 
0.0003 ± 0.0001 
0.00005 ± 0.00002 



TABLE I: Positions of the interfaces and measured values of 
{pi} and {Ui} for the Maier-Stein problem. 

Using FFS, we obtained P B = [4.501 ± 0.007] x 10~ 3 . 
The values of {pi} were also measured (using FFS) and 
are given in Table [I] The landscape variance {Ui} was 
measured using the procedure described in Appendix iDl 
after generating a correctly weighted collection of points 
at interface Xi (for example using FFS), one fires ki trials 
from each point j and records the number of successes, 
Ns \j. One then calculates the variance among points 
V[Nt , ]. The intrinsic variance is given by 



V[N®]/ki-piqi 



(56) 




FIG. 7: Costs of trial runs between interfaces, for the Maier- 
Stein system. The average length, in simulation steps, of 
"successful" trials (to Xi+i) are shown as filled circles. For 
these trials, Xj — Xi + 0.2 and \Xi — Xj\ = 0.2. The average 
length of "unsuccessful" trials (to \a = —0.7) are shown as 
open circles. For these trials, \Xi — Aj| = Xi + 0.7. The solid 
lines show the linear approximation, Eq.JSJ, with S = 131. 



(a maximum of 0.27 for interface 0), indicating that the 
landscape variance is unlikely to have important effects 
in this case. However, this may not be the case for more 
complex systems in higher dimensions. 



B. Testing the expressions 

We now measure directly the cost, in simulation steps, 
the error in the calculated rate constant, and thus the effi- 
ciency of the three methods, for the Maier-Stein problem, 
and compare our simulation results to the predictions of 
Section II I II For each method, simulations were carried 
out in a series of blocks. For FFS, a block consists of a 
complete FFS calculation with Nq starting points. For 
the RB and BG methods, a block consists of No starting 
points at Ao. Each block produces a result P B for the 
probability of reaching B. To find V[Pg], we calculate 
the variance between blocks: 



V[P%] = {PI 



(P L B 



(57) 



Table J] shows that for this problem Ui/pf is rather small 



where the over-line denotes an average over the blocks. 
The cost C per starting point at Ao is the average number 
of simulation steps per block, divided by N$. 

Figure |H1 shows a comparison between the simulation 
values of C and the theoretical predictions (Eqs 10 , i|ll|) 
and H12|) ). for the three methods, as a functions of k 
(FiglHt) and of n (FiglHb). In these calculations, the 
same value of k was used for all interfaces: ki = k for 
all i. To obtain the data in FigHJ), we used interfaces 
which were evenly spaced in A and a fixed value k = 
3. We observe remarkably good agreement between the 
predicted and observed values for the cost, verifying that 
at least for this problem, Eqs 0, l|ll|) and (|12f> are very 
accurate. 

The predicted and measured values of V are shown in 
Figure El for all three methods. Agreement is again ex- 
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FIG. 8: Predicted and measured values of C, for the Maier- 
Stein problem as described in Section II VI The lines show 
the theoretical predictions for the FFS (solid line), BG (dot- 
ted line) and RB (dashed line) methods. The symbols 
show the simulation results. Circles: FFS method, squares: 
BG method, triangles: RB method (with Metropolis accep- 
tance/rejection). Simulation results were obtained with 400 
blocks of No = 1000 starting points for FFS and 2000 starting 
points per block for BG and RB. (a): C as a function of k, for 
n — 7. (b): C as a function of n, for k = 3, for evenly spaced 
interfaces. 
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FIG. 10: Predicted and measured values of V, for the Maier- 
Stein problem, for the FFS method. Solid line: Eg. 1501 (with 
landscape variance), dotted line: Eg. 1231 (no landscape vari- 
ance), circles: simulation results. 
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FIG. 9: Predicted and measured values of V, for the Maier- 
Stein problem. The lines show the theoretical predictions for 
the FFS (solid line), BG (dotted line) and RB (dashed line) 
methods. The symbols show the simulation results. Circles: 
FFS method, squares: BG method, triangles: RB method 
(with Metropolis acceptance/rejection). Simulation results 
were obtained with 400 blocks of No = 1000 starting points 
for FFS and 2000 starting points per block for BG and RB. 
Interfaces were evenly spaced between A a = —0.7 and A_g = 
0.7 (a): V as a function of k, for n — 7. (b): V as a function 
of n, for k — 3. In (b), the landscape contribution is not 
included in the theoretical calculation. 



cellent, showing that the approximations of Section lill Bl 
are justified, at least for this problem. The landscape 
contribution to V is included in Figure [5] for panel (a) 
but not for (b). In Figure ^3 we show the effect of ne- 
glecting this contribution (note the altered scales on both 
axes). Although the landscape contribution is small for 
this problem, it becomes significant for large k as the 
"binomial" contribution decreases. 

The efficiency £ is plotted in Figure ^2 Excellent 
agreement is obtained between simulation and theory. It 
is also interesting to note that the trends in £ as a func- 
tion of k are qualitatively very similar to those obtained 
for the model problem of Fig. [3J The BG method shows 
high efficiency only within a relatively narrow range of 



6 7 
n 



FIG. 11: Predicted and measured efficiency £, for the Maier- 
Stein system. The lines show the theoretical predictions for 
the FFS (solid line), BG (dotted line) and RB (dashed line) 
methods. The symbols show the simulation results. Circles: 
FFS method, squares: BG method, triangles: RB method 
(with Metropolis acceptance/rejection). Simulation results 
were obtained with 400 blocks. For FFS, each block had No = 
1000 starting points and for BG and RB each blocks had 2000 
starting points. Interfaces were evenly spaced, (a): £ vs k for 
n — 7. (b): £ vs n for k = 3. 



parameter values, while the FFS and RB methods are 
much more robust to changes in the parameters. The RB 
method is consistently less efficient than FFS, due to the 
acceptance/rejection step. As the number of interfaces n 
becomes large, we would expect the correlations between 
interfaces (which are not included in our analysis) to have 
a greater effect, and the theoretical predictions to become 
less accurate. This effect is observed to a certain extent: 
the efficiency of FFS, for example, decreases relative to 
the predicted value as n increases. However, this is not a 
dramatic effect, and in fact, even on increasing n further, 
as far as 100 interfaces, we find a decrease of only a few 
percent in the efficiency of FFS. It seems therefore, that 
for FFS at least, one can use any number n of interfaces, 
as long as n is not too small or so very large that memory 
requirements become the limiting factor. 

The remarkable agreement between the theoretical pre- 
dictions and the simulation results shown in Figures |3 
and II II perhaps reflects the simplicity of the Maier-Stein 
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problem. The main assumption for the calculation of V 
- that the sampling of pi at different interfaces is un- 
correlated - seems to be well justified in this case. We 
would expect our theoretical predictions to be less ac- 
curate for more complex problems, perhaps with strong 
correlations between interfaces. In fact, on investigating 
the two examples presented in our previous paper [jj - 
the flipping of a genetic switch and the translocation of 
a polymer through a pore - we find that the quantitative 
estimates of both the cost and variance can differ by a 
factor of about 10 from the theoretical predictions. Even 
with this caveat, however, we believe that the expressions 
of Section ITTll will prove to be of practical use for a wide 
range of rare event simulation problems. 



V. DISCUSSION 

In this paper, we have derived simple analytical expres- 
sions for the computational cost of the three FFS-type 
rare event simulation methods and the statistical accu- 
racy of the resulting estimate of the rate constant. The 
expressions were found to be in remarkably good agree- 
ment with simulation results for the two-dimensional 
non-equilibrium rare event problem proposed by Maier 
and Stein [HEEll- 

Our analysis allows us to draw some general conclu- 
sions about the relative merits of the three FFS-type 
methods. Firstly, the optimum efficiencies of the meth- 
ods are all of the same order of magnitude, at least for 
the simple test problem studied here. However, the meth- 
ods show very different sensitivities to the choice of pa- 
rameters. The Branched Growth method in particular is 
highly sensitive, performing well only for a narrow range 
of parameter values. Within this range, however, it per- 
forms well in comparison to the other methods. The FFS 
method is the most robust to changes in the parameters, 
performing consistently well, even for parameter values 
where the other methods are very inefficient. The Rosen- 
bluth method is lower in efficiency than the others, as a 
consequence of the Metropolis acceptance/rejection step 
which is required in order to obtain paths with the cor- 
rect weights in the Transition Path Ensemble. 

These observations provide a very useful guide for 
choosing a rate constant calculation method. In gen- 
eral, unless one has a very good idea of the optimum 
parameters, the BG method carries a risk of being low 
in efficiency. Of course, strategies could be envisaged to 
overcome this problem - for example, one could imagine 
terminating a certain percentage of the branches to avoid 
the high cost of sampling later interfaces. The analysis 
used here could easily be extended to predict the likely 
success of such approaches. The RB method appears 
from this analysis to be of relatively low efficiency. How- 
ever, that is not to say that one should not use the Rosen- 
bluth method. On the contrary, this is the only method 
which generates unbranched paths, making it highly suit- 
able for situations where one wishes to analyse the paths, 



in order to study the transition mechanism. The RB and 
BG methods also require much less storage of system con- 
figurations than FFS (for which all points at interface 
i must be stored in memory), making them potentially 
suitable for large systems. As a general conclusion, how- 
ever, the results of this paper show that the FFS method 
is highly robust to parameter changes and is probably 
the method of choice for calculations of the rate constant 
where effects such as the storage of many configurations 
in memory are not important. 

These results could also suggest possible strategies for 
choosing the parameters for the three methods. One ap- 
proach would be to use the analytical expressions derived 
here in an optimization scheme for finding {fci}, {A^} and 
n. This is likely to be useful for the BG method, but may 
be less essential for the FFS and RB methods, where the 
choice of parameters is much less critical. 

We expect that the predictions of the cost and statisti- 
cal error derived here will be useful not only for parame- 
ter optimization, but also for assessing, before beginning 
a calculation, which method to use and, indeed, whether 
to proceed at all. Some preliminary calculation would 
be needed in order to obtain rough estimates for 
Pb, {Pi} and (if required) {Ui}. These preliminary cal- 
culations are expected to be much cheaper than a full 
simulation. While the expressions for the cost and vari- 
ance will be less accurate if only rough estimates for the 
parameters are available, we expect the results to be nev- 
ertheless accurate enough to be of use. 

Furthermore, the expressions for V can be used, after a 
rate constant calculation has been completed, to obtain 
error bars on the calculated value of ftAB- In this case, 
the values of Pb and {pi} are known. The intrinsic vari- 
ances {Ui} can also be easily obtained during the rate 
constant calculation, as explained in Appendix^] These 
values can be substituted into the expressions to obtain 
a reliable estimate of the statistical error in the resulting 
rate constant. 

In this work, we provide a way to compare the effi- 
ciency of the three FFS-type methods. It would also be 
very useful to compare their efficiency to that of other 
methods, such as the method of Crooks and Chandler £| 
for non-equilibrium rare event problems, or TPS or 
Transition Interface Sampling (TIS) [a, |2(| for equilib- 
rium problems. We have carried out preliminary calcula- 
tions using the Crooks-Chandler method for the Maier- 
Stein system. We find that the value of the rate con- 
stant is in agreement with that of the FFS-type methods, 
but that the FFS-type methods are much more efficient. 
However, a thorough comparison would require a detailed 
investigation, optimizing the parameter choices of all the 
methods. We therefore leave this to a future study. 

In conclusion, we have presented expressions for the 
computational cost and statistical accuracy of three re- 
cently introduced rare event simulation methods. We 
believe that the expressions presented here will be valu- 
able in using these methods to compute rate constants 
and in evaluating the results of such computations. 
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APPENDIX A: COST OF TRIAL RUNS 



In this case, the average time for a particle to be cap- 
tured at a specified boundary is linearly proportional to 
the distance between the starting point of the particle 
and that boundary, and the proportionality constant is 
the same for particles moving against or with the drift 
velocity. It is therefore appropriate to approximate the 
cost of a trial run between Ai and Xj by S\Xj — Aj|, as in 
Eq.©. 

APPENDIX B: ACCEPTANCE PROBABILITY 
FOR THE RB METHOD 



In order to estimate the cost of a trial run, we as- 
sume that the system undergoes one-dimensional diffu- 
sion along the A coordinate, with a constant drift velocity 
(the origin of which is a force due to the "free energy bar- 
rier"). The problem is then equivalent to that of a parti- 
cle which undergoes diffusion with drift along the x axis, 
after being released between two absorbing boundaries. 
We are interested in the mean time t«- or r_> that the 
particle takes to be captured at the left or right bound- 
ary, given that it is eventually captured at that particular 
boundary. Farkas and Fiilop have studied the problem of 
one dimensional diffusion with drift |28j| . They give ana- 
lytical expressions for the probabilities and that 
the particle is absorbed at the left and right boundaries, 
respectively, and the rates of absorption, j«_ and at 
the left and right boundaries. The mean first passage 
time r is the average time before the particle is absorbed 
at one of the boundaries: 



t [j_ + j^] dt 



(Al) 



To compute t«_ and r_^, we require integrals similar to 
Ea. (|AlJ) . but including only events where the particle 
reaches the desired boundary. The integrals must also be 
normalized by the probability of reaching that boundary: 



71- 



(A2) 



Carrying out the integrals 
Farkas and Fiilop for j^, 
their paper |2j|), we arrive at 



2j using the expressions of 
and (Eqs (3-5) of 
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(1 — a)Lv 
2D 



t coth 



aLv 
~2~D 



(A3) 



where v is the drift velocity, D is the diffusion constant, 
the absorbing boundaries are at x = and x = L and 
the particle is released at x — aL at time t. In the limit 
that the drift velocity is large, cosh [Lv/ (2D)] — > 1 and 
tv_ and t_> reduce to: 



T< — 



olL 

v 



(l-a)L 



(A4) 



This section is concerned with the Metropolis accep- 
tance/rejection step in the Rosenbluth method. We de- 
rive the approximate expression l|44|l for the probability 

9i that a newly generated estimate — iVj , /k% for 

the probability pi is accepted. Upon reaching interface 
i, we calculate the Rosenbluth factor = Ilj=o 

corresponding to the newly generated path leading to 
interface i. We compare this to the Rosenbluth factor 
Wj corresponding to the previous path to have been 
accepted at interface i. Acceptance occurs if the ra- 
tio Zi = /W^ is greater than a random number 
< s < 1. If we know the distribution function P(Zi), 
the acceptance probability is given by: 

/■l /"OO 

9i = ds dZi P{Zi) (Bl) 

JO is 



We would therefore like to calculate P{Zi) = 
P(W^ /W^). To obtain this, we require the distribu- 
tion functions for both and Wl°\ We begin with 
Ws . which we can write as 



log [W { 



(n)l 



j=0 



(B2) 



We assume that the log [A ^ i• : '' , ] for each interface j are in- 
dependent variables (i.e. that the sampling at different 
interfaces is uncorrelated) . Since we are adding many 
independent variables, we apply the Central Limit The- 
orem |2|j to Ea. l|B2|l . In the limit of a large number of 



interfaces, the distribution of y. 



(n) 



'2-k 



exp 



= log [W^ 

fa (n) -^) 2 * 
2a? 



where 



and 



Mi = ^£[logiV»] 

3=0 



i-1 

a\ = ^F[logAT«] 

3=0 



(B3) 



(B4) 



(B5) 
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The expectation value E [log Ns] can be found approx- 

( 7 1 

imately by performing a Taylor expansion of log Ns 
about £[7V S W ], to give: 



(a^-£[a^]) 

logiV« » lo gj B[iVW] + ^ >- (B6) 



( 



N^ - E\N ( S J) 



E[N 

2 



2 £[ivj j) ] 2 
taking the expectation value of Eq.{B6j, we obtain: 



E[logN^] « log£;[JVp)]- (B7) 
Using the variance relation (|17|) . we find that 



F[logiV«] « _l_F[iV«] 



(B8) 



We now need to know E[M j) ] and V[M j) }. On firing 
ki trials from interface i, we know that the number of 
successes follows a binomial distribution. However, the 

variable Ns in Eqs 1|B13(I and l|B14(l refers to the number 
of successes at interface j, given that we know the path 
subsequently reached interface i > j. We therefore know 



that ivi j) > 0, so that 
p{N U)) = _J 



k,\ 



N^>"> ki-N ( J1 
Pi 1, 



(l-q^ikj-N^y.iN^Y' 3 



so that 



25(JV«>) = kjP { 
(l-#) 



(B9) 
(BIO) 

(BH) 



and 



F[iV«] = 



: i % -i^i'i'h '•■ji'j'i'j 



i Y2 



(B12) 



Substituting (|BT0|) and l|B12)l into JE3 and CEHl, w e 

obtain: 



£[logiV«] » log 
and 

V[logN^] 



kjPj 


1 






~ 2 


kjPj 



kjpj 



#3) 



(B14) 



Substituting (|B13|) and (|B14|) in turn into (|B4|) and JB5JI 

leads to 



i-i 

3=0 

and 



k^ 


1 






to | 


kjPj 



1i 



j=0 



(B15) 



(B16) 



Finally, the distribution function f(Wi) for the Rosen- 
bluth factor of the newly generated path can be found 
by making the change of variables W% = exp [yj ] in 
Ea.l|B3j). to give: 



f(Wi) = 



1 



1 



IV, 



exp 



(log[Wi]-/ii) 



2 1 



2a? 



(B17) 



We now turn to the distribution function g(Wi) for the 
Rosenbluth factor of the previous path to have been 
accepted at interface i. does not follow the same 

distribution as W- , because the "previous" path has 
survived at least one round of acceptance/rejection. We 
know that the acceptance/rejection procedure re- weights 
paths by a factor proportional to the Rosenbluth factor 
(see Section HlC(l . so if we assume that has been 

"fully" re- weighted (note that this is an approximation), 
we can say that 



g(Wi) 



f™W'f(W')dW' 



(B18) 



The denominator of Ea. l|B18| > ensures that giWi) is prop- 
erly normalized. Substituting (|B18|I into (|B17|I . we find 
that: 



■ exp 



Oo g rwi]- W 



where 



r-oo 

I = / Wif{Wi)dWi = exp 
Jo 



2a? 



(B19) 



(B20) 



Armed with Eqs (|B17fl and IIB19fl . we can now find 
the distribution function P(Zi) for the ratio Zi = 

Wf/W^K This is given by: 

r°° r°° fw' \ 

P(Zi) = J jf dWidWlgiWi) f{Wl) 6 I ^ - ZA 

(B21) 

Changing the variable of the second integral to Z[ — 
W-/Wi, we obtain 

/>oo poo 

P{Zi) = / / dW l dZ' l W l g{W i )f{Z' i W i )5{Z' l -Z, i ) 
Jo Jo 

dW l W l g{W l )f{Z l W l ) (B22) 
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Substituting 1|B17|) and l|B19jl into l|B22|) . we obtain 
1 



P{Zi) = 



dWi exp 



(B23) 
(log [Z,Wi]- W ) 2 



2o 



This integral can be carried out analytically |30j , to give: 

ex P^-^l r (io g z 4 ) 2 log Zi 



P{Zi) = 



2a iZi^/n 



exp 



4(7? 



(B24) 



We are now finally in a position to calculate the ac- 
ceptance probability 9 i: using Ea. (|Bl|) . Substituting 
Ea. (IB24ll into (|B1|) and integrating over Zi, we obtain 
EH: 





1-^Erf 


Pi 


logs" 




f ds 








la 
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2 " 


2a t _ 





(B25) 



2Erf 



(1± 
\ 2 



- 1 



Erf(x) 



where Erf (2;) is the error function: 
(2/V^)Joe- t2 dt. 

Although Eo. l)B25|l is a simple and convenient expres- 
sion for the acceptance probability 9i, its derivation re- 
quired several approximations. We have therefore tested 
the validity of Ea. (|B25|l . We first carried out a "simu- 
lated simulation" , in which we defined a series of N = 15 
interfaces, each with the same value of pi = p = 10~ 6 / 15 , 
and "simulated" the Rosenbluth calculation, each time 
drawing a random number to determine the outcome of a 
given "trial run" , for a given number of trial runs fcj = k, 
taken to be the same for all interfaces. We measured the 
acceptance probabilities at each interface after 2 x 10 6 
Rosenbluth "path generations" , and compared these to 
Ea. (IB25|l . The results are shown in FigurefHlk. for k = 2, 
k = 5 and k = 8. The agreement with the "simulation" 
is very reasonable. To compare with real simulation re- 
sults, we also measured the acceptance probabilities 9i for 
the RB simulations of the Maier-Stein system described 
in Section llVl The results are compared with the pre- 
dictions of Ea. (|B25|) in Figure fElb. Again, quite good 
agreement is obtained. 



APPENDIX C: THE EFFECTS OF LANDSCAPE 
VARIANCE 

In this section, we include the effects of the "landscape 
variance" in our expressions for the relative variance V 
of Pg. The result will be that expressions (|23|l . 
and |g3J are transformed into (J5TJ and As 

described in Section IIIII we suppose that point j at in- 

(7) 

tcrface Ai has probability pf that a trial run fired from 
it will reach Aj+i, rather than Xa- The variance in the 
Pi values for points at Ai (sampled according to their 
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FIG. 12: (a): "Simulated" and predicted acceptance proba- 
bilities 8i for interfaces < i < 14, for the "simulated simula- 
tion" described in the text, for k — 2, 5, 8. (b): Simulated and 
predicted values of Oi for < i < 6 for the Maier-Stein prob- 
lem of Section II VI for k = 2,5, 8. In both plots, solid lines 
represent predicted values for k — 2, dotted lines, k = 5 and 
dashed lines, k — 8. Symbols represent simulation results: 
circles: k — 2, squares: k — 5 and triangles: k — 8. 



expected occurrence in the trial run firing procedure) is 
the "landscape variance" , C/j. 

If we choose a particular point j, fire ki trial runs 



and measure the number of successes N. 



to obtain a mean value E[N l f , \j] = hp 



we expect 
and a vari- 



ance V[iVi^|j] = fcipp 'q'f' (in analogy with Eqs (|T5|) and 
p6[lh We now average over many points j at interface 
Aj, using the general variance relation 12(il) : 



V[N&] = E V N^\j 



V 



E 



NP\j (CI) 



E Ihp^q^} + V \hp 



where the mean and the variance are taken over the 



distribution of points j. Since E\p^q^] = E[p, 



0')\2l _ 

„0O 



E\p\ j) ] - E[{p^) 2 ] and U t = E[{p 



(P.. 

(Elpf'}) 2 , we can deduce that E kipf'q. 



C7')\2l 



= h(pi - 

Ui) = ki(piqi — Ui). Using Ea. itH|) . we have 



V[kipP] = kfV^p] = kpi, so that 
V[N®} = k tPl q t + Ui& 



(C2) 



This first term on the r.h.s. of Ea. (|C2(l corresponds to 
Ea. (|16|) : the binomial contribution arising from the lim- 
ited number of trial runs per point. The second term is 
an extra contribution, due to the landscape variance. 

We now repeat the derivations of Section fill Bl simply 
replacing Eq.JTBJl by Ea. l(U2|) . We begin with the RB 
method, for which Ea. (|41|) becomes 



' 1 " 




~p%q% 


No. 







Ui 1- 



(C3) 



(2-0* 



ni-o('-^) 
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and Eq. {i^ is replaced by Eq. : 



i=0 



(2-0* 



For the BG method, Ea. lp^Hjl is replaced by: 

= [**Pi«i + ^ (*? - **)] nto *iPj 



+kfp 2 i V[Nj S i - 1) ] 
= kiPiqi + Ui (fcf - hi) 

and Eq.J23) becomes Ea. (|5T"]) : 

i=o '' " 1 



(C4) 

(<>0) 
(i = 0) 



v bs = E 



For the FFS method, the situation is slightly more 
complicated, because the number of trials fired from 
point j at interface i is not fixed. We make Mi trials 
from the N points at Aj, each time selecting a start- 
ing point at random (so that the probability a particular 
point is chosen is X/Nf). Since we no longer assume that 
all points at interface i are identical, we must now take 
account of the distribution of the number of times rrij 
that point j is selected. This is in fact a multinomial 
distribution H^,|3^|, which has average E[rrij] = Mi/Ni 
and variance U[mj] = Mi [l/iVj(l — 1/Nf)]. Let us now 
do a "thought experiment" in which we first decide how 
many trial will be fired from each point - i.e. we fix the 
set of values {rrij} (of course, J2j m j = Mi). We then 
fire these trials and measure the total number 7V* ot which 
reach A»+i. The expectation value for Nl ot is 



E m irt 

3 



M iPi 



(C5) 



and the variance is found using Ea. (|C2(l . with k% replaced 
by rrij , multiplying by m 2 and summing over all j : 



V[N t s ot \{m j }] = J2[m j PiQi + Ui[r 



(C6) 



We now imagine that we average the results over many 
sets of values {rrij}. Using the general relation l(2Bjl. we 
obtain: 

V[N?*\ = V[E[Nl ot \{ mj }]] + E[V[Nt ot \{m 3 }]] (C7) 



= V[MiPi]+E 



MiPiqi 



^E m ? 



U,M t 



= Mip iqi + Ui [NiE[mj] - Mi] 

Here, the variance and expectation are with respect to 
the distribution of {rrij} values. The last line follows 



from the fact that V [MiPi] — as both Mi and Pi 
are constants with respect to changes in {rrij}. Since 
V[m 3 ] = M t [l/Ni(l - 1/Ni)] = E[m 2 j] - E[rrij} 2 , we find 
that E[m 2 } = (Mi/N t ){l-l/N t ) + Mf/Nf. Substituting 
this into Ea. (|C7fl . we obtain: 

V[Nt ot ] = Mi^qi + ^ [Mf - Mi]] (C8) 



elffs. 



Since p\ = Nl ot /M u we must divide Eq.(C8j! by Mf to 
obtain V\p\ 



ifis _ Pili 



Ui 



This leads to: 



n-l 



V Ss = N J2 



Mi Ni 



q, 



l 

M 



(C9) 



2 = 



U 



Pi Mi P jN t 

n;= (i-^) 



1 

M 



(CIO) 



where Ni = Mi_ipj_i for i > and N, = Nq for i = 0. 
Rewriting in terms of fe$ = Mi/No, we obtain Ea. (|50f) : 



V 



lis 



E 

i=0 



Qi 


, U t N 


Pih 


+ P 2 N 














nj= 





1 - 



APPENDIX D: MEASURING THE INTRINSIC 
VARIANCE 




FIG. 



13: V[N^ 0) ]/k 2 (solid line) and (l/(fc - 
1)) [u[ivi 0) ]/fc — po<7o (dashed line), as functions of 
k = Mo /No, calculated using FFS as described in Section IHl 
for the Maier-Stein problem of Section llVl with 10000 points 
at the first interface Aq = —0.7. 
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In this section, we describe a simple and computation- 
ally cheap procedure for measuring the landscape vari- 
ance parameters C/j. Given a correctly weighted collec- 
tion of Ni points at interface Xi (obtained, for example, 
using FFS), we could fire an extremely large number k 
of trial runs from each point and measure the variance 
among points in the values of N^'^ - where Ng 1 '^ de- 
notes the number of successful trials from point j: 



U % = V[ P i 



V[N, 



(Oi 



A- 2 



3=1 



N, 



^A 



3=1 



N, 



( D1 ) 

This is likely to be an expensive procedure. Fortunately, 
however, it is not necessary to fire a very large number 
of trial runs from each point. Instead, we make use of 
expression (IC2JI . which can be written as 



Ui = 



kV{Pl 



1 



k - i 



(fc-i) 



(D2) 



where the expression now holds for any value of k. In 
the limit that k — ► oo, Ea. (|D2|l reduces to (|D1|) . As a 
practical procedure, therefore, we generate a collection 
of Ni points at interface i (using, for example, FFS), 
and fire k trials from each point - k does not have to be 
a large number. For each point j, we record the num- 



ber of successful trials Ns ■ The variance v[m i> ] of 



Wl 



these values is inserted into Ea. (|D2|l to give a value for 
Ui. Figure IT51 shows the results of this procedure for the 
Maier-Stein problem of Sect ion HVI For the first interface 
(Ao (S =$&Q), Ui was calculated using Eq.([D2j|, using k 
trials for each of 10000 points collected at Ao- The solid 
line is the measured value of V[Ns ]/k 2 , while the dashed 
line is the value of Ui obtained from Ea. (|D2|) . The two 
lines converge, of course, for large values of k. Figure ITS1 
shows that accurate results for C/j can be obtained us- 
ing Eq. (|D2() , using only a small number of trial runs per 
point. 
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